#include "../include/Matrice.h"
#include <math.h>
#include <stdlib.h>
#include <iostream>

Matrice::Matrice ( int _largeur, int _hauteur )
{
	this->largeur = _largeur;
	this->hauteur = _hauteur;
	this->matrice = new float*[_largeur];
	for ( int i = 0 ; i < _largeur ; i++ )
	{
		this->matrice[i] = new float[_hauteur];
	}
}

Matrice::~Matrice ()
{
	//std::cerr << "destructeur de matrice" << std::endl;
	for ( int i = 0 ; i < largeur ; ++i )
	{
		delete this->matrice[i];
	}
	delete this->matrice;
}

Matrice::Matrice ( int _largeur, int _hauteur, float** _matrice )
{
	this->largeur = _largeur;
	this->hauteur = _hauteur;
	matrice = _matrice;
}

Matrice::Matrice ( int _largeur, int _hauteur, float* _matrice )
{
	this->largeur = _largeur;
	this->hauteur = _hauteur;
	matrice = new float*[_largeur];
	
	for ( int i = 0 ; i < _largeur ; i++ )
		matrice[i] = new float[_hauteur];
	
	for ( int i = 0 ; i < _largeur ; i++ )
		for ( int j = 0 ; j < _hauteur ; j++ )
			this->matrice[i][j] = _matrice[j*_largeur+i];
}

float Matrice::getValeur  ( int i, int j )
{
	if ( i < 0 )
		i = 0;
	if ( i >= largeur )
		i = largeur-1;
	if ( j < 0 )
		j = 0;
	if ( j >= hauteur )
		j = hauteur-1;
	return this->matrice[i][j];
}

void Matrice::setValeur  ( int i, int j, float v )
{
	if ( i < 0 )
		i = 0;
	if ( i >= largeur )
		i = largeur-1;
	if ( j < 0 )
		j = 0;
	if ( j >= hauteur )
		j = hauteur-1;
	matrice[i][j] = v;
}

float Matrice::getLargeur  ( )
{
	return this->largeur;
}

float Matrice::getHauteur  ( )
{
	return this->hauteur;
}


float Matrice::trace()
{
	if ( this->hauteur != this->largeur )
	{
		std::cout<<"t'as essaye nawak avec tes matrices la, je me barre\n";
		exit(1);
	}
	
	float somme=0.0;
	
	for ( int i = 0 ; i < largeur ; i++ )
		somme += matrice[i][i];
	return somme;
			
}

Matrice* Matrice::multiplication ( Matrice* m, Matrice * resultat/*=NULL*/ )
{
	if ( m->hauteur != this->largeur )
	{
		std::cout<<"Et voila encore de la merde\n";
		this->afficher();
		std::cout<<"\n\n";
		m->afficher();
		exit(1);
	}
	float somme;
	
	if ( resultat == NULL )
		resultat = new Matrice (m->largeur,this->hauteur);
	/*std::cout<<"mat resultat > largeur et hauteur : "<<resultat->largeur<<" "<<resultat->hauteur<<std::endl;
	std::cout<<"mat this > largeur et hauteur : "<<this->largeur<<" "<<this->hauteur<<std::endl;
	std::cout<<"mat m > largeur et hauteur : "<<m->largeur<<" "<<m->hauteur<<std::endl;*/
	//std::cerr << "*c" << std::endl;
	/*for ( int i = 0 ; i < m->largeur ; i++ )
	{
		for ( int j = 0 ; j < this->hauteur ; j++ )
		{
			somme = 0;
			std::cerr << "*d" << std::endl;
			for ( int k = 0 ; k < m->largeur ; k++ )
			{
				std::cerr << "*e" << std::endl;
				somme += this->matrice[i][k] * m->matrice[k][j];
				std::cerr << "*f" << std::endl;
			}
			
			resultat->matrice[i][j] = somme;
		}
	}*/
	
	for ( int x = 0 ; x < resultat->largeur ; ++x )
	{
		for ( int y = 0 ; y < resultat->hauteur ; ++y )
		{
			somme = 0.0;
			for ( int xx = 0 ; xx < this->largeur ; ++xx )
			{
				somme += this->matrice[xx][y] * m->matrice[x][xx];
			}
			resultat->matrice[x][y] = somme;
		}
	}
	
	return resultat;
}

Matrice* Matrice::multiplication ( float _x, Matrice * resultat/*=NULL*/ )
{
	if ( resultat == NULL )
		resultat = new Matrice(this->largeur,this->hauteur);
	
	for ( int i = 0 ; i < this->largeur ; i++ )
		for ( int j = 0 ; j < this->hauteur ; j++ )
			resultat->matrice[i][j] = this->matrice[i][j] * _x;
	
	return resultat;
}		

Matrice* Matrice::addition  ( Matrice* m, Matrice * resultat/*=NULL*/ )
{
	if (  this->hauteur != m->hauteur || this->largeur != m->largeur )
	{
		std::cout<<"Erreur lors de l'addition de matrices\n";
		exit(1);
	}
	
	if ( resultat == NULL )
		resultat = new Matrice(m->largeur,m->hauteur);
	
	for ( int i = 0 ; i < m->largeur ; i++ )
		for ( int j = 0 ; j < m->hauteur ; j++ )
			resultat->matrice[i][j] = this->matrice[i][j] + m->matrice[i][j];
	return resultat;
}

Matrice* Matrice::soustraction  ( Matrice* m )
{
	if (  this->hauteur != m->hauteur || this->largeur != m->largeur )
	{
		std::cout<<"Erreur lors de la soustraction de matrices\n";
		exit(1);
	}
	
	Matrice* resultat = new Matrice(m->largeur,m->hauteur);
	
	for ( int i = 0 ; i < m->largeur ; i++ )
		for ( int j = 0 ; j < m->hauteur ; j++ )
			resultat->matrice[i][j] = this->matrice[i][j] - m->matrice[i][j];
	return resultat;
}

float Matrice::valeurPropreMax ()
{
	float B,C;
	float delta;
	float x1,x2;
	
	/* B et C sont les coefficient du polynome caracteristique de la matrice */
	B = -this->matrice[0][0]-this->matrice[1][1];
	C = this->matrice[0][0]*this->matrice[1][1] - this->matrice[0][1]*this->matrice[1][0];
	
	delta = B*B - 4*C; //(A vaut 1 toujours)
	
	if ( delta < 0 )
	{
		std::cout<<"discriminant negatif, Erreur\n";
		exit(1);
	}
	
	x1 = ( -B -sqrt(delta) )/2;
	x2 = ( -B + sqrt(delta) )/2;
	
	if ( x1 > x2 )
		return x1;
	else
		return x2;
}

float Matrice::valeurPropreMin ()
{
	float B,C;
	float delta;
	float x1,x2;
	
	/* B et C sont les coefficient du polynome caracteristique de la matrice */
	B = -this->matrice[0][0]-this->matrice[1][1];
	C = this->matrice[0][0]*this->matrice[1][1] - this->matrice[0][1]*this->matrice[1][0];
	
	delta = B*B - 4*C; //(A vaut 1 toujours)
	
	if ( delta < 0 )
	{
		std::cout<<"discriminant negatif, Erreur\n";
		exit(1);
	}
	
	x1 = ( -B -sqrt(delta) )/(float)2;
	x2 = ( -B + sqrt(delta) )/(float)2;
	
	if ( x1 < x2 )
		return x1;
	else
		return x2;
}

int Matrice::valeursPropres ( float * _lambda1, float * _lambda2 )
{ // si delta == 0, on aura *_lambda1 == *_lambda2
	if ( largeur != 2 || hauteur != 2 )
	{ // ne fonctionne que pour les matrices 2x2
		std::cout << "Matrice::valeursPropres () : erreur : ne fonctionne que pour les matrices 2x2" << std::endl;
		return -1;
	}
	
	if ( _lambda1 == NULL || _lambda2 == NULL )
	{
		std::cout << "Matrice::valeursPropres () : erreur : pointeur nul donné en paramètre" << std::endl;
		return -1;
	}
	
	float B,C;
	float delta;
	
	/* B et C sont les coefficient du polynome caracteristique de la matrice */
	//B = -this->matrice[0][0]-this->matrice[1][1]; // -a-d
	B = this->matrice[0][0]+this->matrice[1][1]; // a+d
	C = this->matrice[0][0]*this->matrice[1][1] - this->matrice[0][1]*this->matrice[1][0]; // a*d - bc
	
	delta = B*B - 4*C; //(A vaut 1 toujours)
	
	if ( delta < 0 )
	{
		std::cout << "Matrice::valeursPropres () : erreur : discriminant negatif" << std::endl;
		return -1;
	}
	delta = sqrt ( delta );
	*_lambda1 = 0.5*( B - delta );
	*_lambda2 = 0.5*( B + delta );
	return 0;
}


Matrice* Matrice::transposee ( )
{
	Matrice* resultat = new Matrice(this->hauteur,this->largeur);
	
	for ( int i = 0 ; i < this->hauteur ; i++ )
		for ( int j = 0 ; j < this->largeur ; j++ )
			resultat->matrice[i][j] = this->matrice[j][i];
	
	return resultat;
}

Matrice* Matrice::transposee ( Matrice* m )
{
	Matrice* resultat = new Matrice(m->hauteur,m->largeur);
	
	for ( int i = 0 ; i < m->hauteur ; i++ )
		for ( int j = 0 ; j < m->largeur ; j++ )
			resultat->matrice[i][j] = m->matrice[j][i];
	
	return resultat;
}

float Matrice::determinant()
{	
	if ( this->largeur != this->hauteur || (this->largeur != 2 && this->largeur != 3) )
	{
		std::cout<<"Erreur on essaie de calculer le determinant d'une matrice pas carée 3x3\n";
		exit(1);
	}
	
	if ( this->largeur == 2 )
	{
		return ( this->matrice[0][0]*this->matrice[1][1] - this->matrice[1][0]*this->matrice[0][1] );
	}
	
	if ( this->largeur == 3 )
	{
		return this->matrice[0][0]*this->matrice[1][1]*this->matrice[2][2] + this->matrice[0][1]*this->matrice[1][2]*this->matrice[2][0] + this->matrice[0][2]*this->matrice[1][0]*this->matrice[2][1]
			- this->matrice[2][0]*this->matrice[1][1]*this->matrice[0][2] - this->matrice[2][1]*this->matrice[1][2]*this->matrice[0][0] - this->matrice[2][2]*this->matrice[1][0]*this->matrice[0][1];
	}
	
	return 666;
}

Matrice* Matrice::vecteurPropreMax()
{
	float lambda;
	Matrice* resultat = new Matrice(1,2);
	
	lambda = this->valeurPropreMax();
	
	resultat->matrice[0][0] = this->matrice[0][0]-lambda;
	resultat->matrice[0][1] = this->matrice[1][0];
	
	return resultat;
}

Matrice* Matrice::vecteurPropreMin()
{
	float lambda;
	Matrice* resultat = new Matrice(1,2);
	
	lambda = this->valeurPropreMin();
	
	resultat->matrice[0][0] = this->matrice[0][0]-lambda;
	resultat->matrice[0][1] = this->matrice[1][0];
	
	return resultat;
}

Matrice * Matrice::vecteurPropre ( float _lambda, Matrice * _res/*= NULL*/,
	int * _err/*= NULL*/ )
{
	if ( largeur != 2 || hauteur != 2 )
	{ // ne fonctionne que pour les matrices 2x2
		std::cout << "Matrice::vecteurPropre () : erreur : ne fonctionne que pour les matrices 2x2" << std::endl;
		return _res;
	}
	if ( _res == NULL )
		_res = new Matrice ( 1, 2 );
	
	float a = getValeur ( 0, 0 ) - _lambda;
	float b = getValeur ( 1, 0 );
	float c = getValeur ( 0, 1 );
	float d = getValeur ( 1, 1 ) - _lambda;
	float precision = 0.00001;
	
	if ( _err != NULL )
			*_err = 0;
	
	if ( (a > precision || a < -precision) && (b > precision || b < -precision) )
	{
		_res->setValeur ( 0, 0, 1.0 );
		_res->setValeur ( 0, 1, -a/b );
	}
	else if ( (c > precision || c < -precision) && (d > precision || d < -precision) )
	{
		_res->setValeur ( 0, 0, 1.0 );
		_res->setValeur ( 0, 1, -c/d );
	}
	else
	{
		//std::cout << "Matrice::vecteurPropre () : attention : valeurs proches de 0" << std::endl;
		_res->setValeur ( 0, 0, 0.0 );
		_res->setValeur ( 0, 1, 0.0 );
		if ( _err != NULL )
			*_err = -1;
	}
	
	return _res;
}

int Matrice::calculValEtVecPropres ( float * _lambda1, float * _lambda2,
	Matrice * _theta1, Matrice * _theta2 )
{
	int res = 0;
	if ( largeur != 2 || hauteur != 2 )
	{ // ne fonctionne que pour les matrices 2x2
		std::cout << "Matrice::calculValEtVecPropres () : erreur : ne fonctionne que pour les matrices 2x2" << std::endl;
		return -1;
	}
	if ( _lambda1 == NULL || _lambda2 == NULL || _theta1 == NULL || _theta2 == NULL )
	{
		std::cout << "Matrice::calculValEtVecPropres () : erreur : pointeur nul donné en paramètre" << std::endl;
		return -1;
	}
	res = valeursPropres ( (float*)_lambda1, (float*)_lambda2 );
	if ( res != 0 )
		return -1;
	
	/*vecteurPropre ( *_lambda1, (Matrice*)_theta1, &res );
	if ( res != 0 )
		return -1;
	vecteurPropre ( *_lambda2, (Matrice*)_theta2, &res );
	if ( res != 0 )
		return -1;
	remplacé par : version de CImg::eigen ()*/
	float t = *_lambda1;
	*_lambda1 = *_lambda2;
	*_lambda2 = t;
	const double theta1 = atan2( *_lambda2 - getValeur ( 0, 0 ), getValeur ( 0, 1 ) ),
		theta2 = atan2 ( *_lambda1 - getValeur ( 0, 0 ), getValeur ( 0, 1 ) );
	_theta1->setValeur ( 0, 0, cos ( theta1 ) );
	_theta1->setValeur ( 0, 1, sin ( theta1 ) );
	_theta2->setValeur ( 0, 0, cos ( theta2 ) );
	_theta2->setValeur ( 0, 1, sin ( theta2 ) );
	
	
	return 0;
}

void Matrice::normaliserVecteur ()
{
	float norme = 0.0;
	for ( int i = 0; i < hauteur; ++i )
	{
		norme += getValeur ( 0, i )*getValeur ( 0, i );
	}
	norme = sqrt ( norme );
	for ( int i = 0; i < hauteur; ++i )
	{
		setValeur ( 0, i, getValeur ( 0, i ) / norme );
	}
}

void Matrice::afficher()
{
	for ( int i = 0 ; i < this->hauteur ; i++ )
	{
		for ( int j = 0 ; j < this->largeur ; j++ )
		{
			std::cerr<<this->matrice[j][i] << "  ";
		}
		
		std::cerr<<std::endl;
	}
}

Matrice * Matrice::mmT ( Matrice * _res/*= NULL*/ )
{
	Matrice * mT = transposee ();
	_res = multiplication ( mT, _res ); // si _res était NULL, il devient le pteur vers la new mat
	delete mT;
	return _res;
}

Matrice * Matrice::mmT2 ( Matrice * _res/*= NULL*/ )
{
	if ( /*largeur != 1 ||*/ hauteur != 2 )
	{
		std::cerr << "Matrice::mmT2 () : erreur : mauvais format de matrice en entrée" << std::endl;
		return _res;
	}
	if ( _res == NULL )
		_res = new Matrice ( 2, 2 );
	else if ( _res->largeur != 2 || _res->hauteur != 2 )
	{
		std::cerr << "Matrice::mmT2 () : erreur : mauvais format de matrice en entrée" << std::endl;
		return _res;
	}
	float a = getValeur ( 0, 0 );
	float b = getValeur ( 0, 1 );
	_res->setValeur ( 0, 0, a*a );
	_res->setValeur ( 1, 0, a*b );
	_res->setValeur ( 0, 1, b*a );
	_res->setValeur ( 1, 1, b*b );
	return _res;
}
